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INTRODUCTION 


There  are  many  instances  in  operations  research  and  statistics 


in  which  one  must  integrate  a  bivariate  normal  density  over  an  offset 


circle.  Namely,  let 


f(x,y)  = 


2  7TCT  a  (1-p  ) 
x  y 


exp 


2  ( 1-p  ) 


(x-m  ) (y-m  ) 
x  J  y 

a  a 
x  y 


(1) 


be  the  density  for  a  bivariate  random  variable  (X,Y) .  The  probability 
PM  ,  that  (X , Y)  falls  in  the  offset  circle  R  defined  by: 


(2) 

(3) 

Equation  3  is  the  coverage  integral  or  coverage  function  (e.g., 
see  references  1  through  6) . 

The  coverage  integral  has  been  calculated  in  certain 
special  cases  (references  6  and  7).  When  m  =  m  =0  and  the 

X  y 

circle  is  centered  at  the  origin,  one  can  convert  the  integrand 


2  2  2 

R  =  ( x , y :  (x-xc)  +  (y-yc)  s  r  } 


is 


P 


Rj  =  /  / f  (x,y)dxdy  . 


to  a  Bessel  function,  expand  in  a  power  series  and  integrate  term 
by  term  (reference  6) . 

If  a  =  a  =  a  ,  then  Pl_rJ  can  be  given  in  terms  of 

X  y  3 

2 

the  non-central  x  distribution  (reference  7) . 

The  result  of  reference  6,  however,  converges  very  slowly 

if  r //la,  >  1  ,  where  a.  =  min (a  ,a  ).  The  result  of  reference  7 
1  1  x  y 

is  not  valid  for  non-circular  densities. 

4 

in  this  paper,  a  formula  is  derived  for  the  calculation 
of  the  coverage  integral  in  the  most  general  case.  A  sample  problem 
is  presented  and  in  the  appendix  a  Fortran  code  is  listed  for 
implementing  the  formula  on  a  computer. 


THE  ITERATED  INTEGRAL 


With  reference  to  figure  1,  it  is  seen  that  equation  3  can 
be  rewritten  as: 


xc+r  yc+^2-(x-xc)2)h 


PCRJ 


-/  f 


.  ,  2  ,  .  2.  *5  Hx»y)  dydx 

xc-r  yc-(r  -(x-xc)  ) 


(4) 


2  2  h 

X  +r  y  + (r  -(x-x  )  ) 

C  jm  C  W 


5  / 


2lTax0y  (1_p 


-  5  J  J  .  2  .  >2.h 

xc_r  yc_(r  ■(x_xc)  > 


exp 


2  ( 1—  p  2 ) 


(=?)’  •  “Rf?)-  R 


jdydx 


(5) 


Let  I (R)  denote  the  integral  in  equation  5,  so  that: 


P^Rj  = 


2Troxay(1_p  ] 


I  (R) 


This  integral  can  be  rewritten  as: 

HR)  =  / 


x  +r 
c 


exp 


xc'r 


exn 


1  /X~mx\ 
2  ( 1-p2)  \  °x  / 


2  2  h 

y  +(r  -(x-x  ) 

f 

yc-(r2-(x-xc)2)h 


ix|Rf  ’XRX 


y-m 

: _ i 

a 


2  (1-p  ) 


dy 


dx 


(6) 


(7) 
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From  equation  11,  one  sees  that: 
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In  terms  of  this  function,  equation  15  is  rewritten  as: 

1,W  -  {-r  4  ift) 

jay  [2  ( 1~p2 )]  Erf  (s+(x,r,xc,yc))-Erf  (s^(x,r,xc,y c)) 


Now  introduce  a  new  variable,  v  ,  so  that 


i.e. ,  dx  =  <Txdv 


and  x  =  vo  +  m 

x  x 


When  x  =  x  +r  ,  v  =  (x  +r~m  )/a  and  when  x  =  x  -r , 

C  c  X  X  c 

v  =  (x  -r-m  )/a  .  Equation  17  can  now  be  rewritten  as: 

^  XX 


x  +r-m 
c  x 


KRI  -  [2(l-p2)]  <Jxo 


2.1'  „  /I  (  °x 

>J  °x°y  ""2  ./  _  _ 


-I''2 


x  -r-m 
c  x 


j^Erf  (s+(axv+mx,r,xc,yc))  -Erf(s_(oxv4-mx,r,xc,yc))j  dv. 


Equation  19  is  still  exact. 
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APPROXIMATION  OF  THE  ERROR  FUNCTION 


In  order  to  simplify  equation  19,  an  approximation  for  the 
error  function  will  be  used.  From  reference  7,  page  229: 


Set  t  =  l/(l+px)  where  p  =  .3275911.  Then 


2  5 

Erf  (x)  =  1  -  e  X  *  2  a.t1  +  e (x) 

i=l  1 


In  this  equation 


a1  =  .254829592 
a2  =  -.284496736 
a3  =  1.421413741 


a4  =  -1.453152027 

ac  =  1.061405429 
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The  error  term  e  (x)  is  such  that: 


c (x)  £1.5x10 


Now  set 


t  +  ( v )  = 


1  +  pSt (cxv+mx'r,xc'yc) 


Using  equation  23  in  equation  19  and  then  equation  19  in  equation  16, 


leads  to  the  following  result. 


8 


PM 


x  +r-m 
c  x 


C  X 


a 

x 


5 

a .  t  (v)1 

i  =  l 


-S 


-  e 


5 

\~> 

i—J 

i=l 


a^t+ (v) 


dv  +  e. 


(24) 


The  error  term  is  bounded  as  follows: 


::  +r-m 
c  x 


s  —  (3  x  10"7) 
/2tt 


(25) 


(26) 


A  finer  bound  on  is: 


£  (3  x  10~7)  \ 


(27) 


( 
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A  SAMPLE  CALCULATION 


In  this 
parameters  were 


section , 
chosen : 


a  sample  problem  is  studied. 


The  following 


ax  =  200 
mx  =  0 


These  parameters  were  chosen  because  the  method  of  reference  6 
can  be  used  to  obtain  p£r_j  in  this  case. 

According  to  reference  6,  set 

a-j_  =  min(ax,ay)  a 2  =  max(ox,ay) 


Then : 


o .  °° 

m  ^  e 


(-1) 


2  2\k  +  1 

2/2o\) 


k=0 


(k+1 )  ! 


f  V>  /  l)1  (k\  (2i\  (,  2.  2V\ 


Equation  30  is  valid  only  for  m  =  m  =  p  =  0  . 

X  X 


(28) 


(29) 


(30) 
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The  calculations  were  done  on  a  Burroughs  6700  computer. 

The  integral  (equation  24)  was  evaluated  by  the  trapezoid  rule. 

This  technique  was  chosen  instead  of  more  advanced  schemes,  because 
it  is  so  simple. 

In  table  1,  results  obtained  using  equations  24  and  30  are 
compared.  The  integration  range  in  equation  24  was  divided  into  150 
steps.  The  accuracy  of  the  two  calculations  is  good,  even  at 
r  =  275.  However,  when  the  series  is  used,  the  amount  of  computer 
time  used  increases  rapidly  with  r  .  For  example,  for  r  =  250, 
it  took  4.2  seconds  to  evaluate  the  series  but  only  0.3  seconds  to 
evaluate  the  integral. 
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TABLE  1 


COMPARISON  OF  METHODS  FOR  EVALUATION 
OF  THE  COVERAGE  INTEGRA!, 


r 

By  Series  (30) 

By  Integral  (24) 

Percent 

Difference 

25 

.030243 

. 030237 

.021 

50 

.110246 

. 110222 

.  022 

75 

. 215952 

.215898 

.  025 

100 

. 325350 

.325261 

.  028 

125 

.426700 

.426570 

.031 

150 

. 516830 

. 516657 

.033 

175 

.596196 

. 595984 

.  036 

200 

. 665824 

.665578 

.037 

225 

.  726484 

.726212 

.037 

250 

. 778770 

.778480 

.  037 

275 


823259 


822962 


036 


REFERENCES 


Guenther,  W.C.  and  P.J.  Terragno,  "A  Review  of  the  Literature 
on  a  Class  of  Coverage  Problems,"  Annals  Math.  Stat.,  1964, 
35:232-260. 

Center  for  Naval  Analyses,  Professional  Paper  221,  "On 
Approximating  the  Circular  Coverage  Function, "Maurice 
Mizrahi,  Unclassified,  Feb  1978. 

Gilliland,  D.C.,  "Integral  of  the  Bivariate  Normal  Distribution 
Over  an  Offset  Circle,"  J.A.S.A. ,  1962,  57:758-768. 

Oberg,  E.N.,  "Approximate  Formulas  for  the  Radii  of  Circles 
Which  Include  a  Specified  Fraction  of  a  Normal  Bivariate 
Distribution,"  Am . Math . Stat . ,  1947,  18:442-447. 

Guenther,  W.C.,  "A  Generalization  of  the  Integral  of  the 
Circular  Coverage  Function,"  Amer.  Math.  Mo.,  1964,  71:278-283. 

Center  for  Naval  Analyses,  OEG  Study  626,  "Probabiltiy  of 
Damage  Problems  of  Frequent  Occurence,"  1959. 

Abramoroitz,  M.  and  I.  Stegun,  "Handbook  of  Mathematical 
Functions,"  NBS ,  Washington,  D.C.  (1964),  QA3  05. 


13 


I 


APPENDIX 


COMPUTER  CODE  IN  THE  GENERAL  CASE 

The  following  Fortran  code  can  be  used  to  evaluate  the 
coverage  integral  in  the  general  case. 


:•  THIS  SUtiFCUIINE  INTEGRA  ITS  i  I  V  A  R I  A  TE  NORMAL  OVER 

C*  AN  OFFSET  CIRCLE.  SX  *  SY  ARE  THE  ST  ANCARO  DEVIATIONS 
C*  OF  THE  CISTRieUI ION  FUNCTION  CPNTPRFD  *T  UX,  UY.  RG 
C*  IS  THE  CCNRELATICN  COEFFICIENT.  THE  INTEGRAL  IS  EVAi_~ 

C*  GATED  OVER  A  CHCLE  OF  RADIUS  $  Cr  N  Tr  RED  AT  X,  Y  USING 
C*  THE  TRAPEZOIDAL  METHOD.  PR  CONTAINS  THE  VALUE  CF  T  Hc 
C*  INTEGRAL. 

S  L  5  R  CUT  INF  3lVlNT(SX»St#UX»LY,RC*K#<*Y.PR) 

VL=(X-R-UX)/SX 

DL=2.*R/SX/2CC. 

A  1=  .  2  54  8?  S5  S  ? 

A2=-.28449E736 
43=1 .4214 1 3741 
A4=-l.453152C27 
A5=l.C614C542S 
P=. 327591 1 

SG=S0RT(1./C?.*(l-R0*RC))) 

PR=C.C 

00  1C  J  =  t  *  2  C  1 

V  =  C J-l  )*OL* VL 

S1=SCRT(R**2-(SX*V*UX-X)«*2) 

3U=3G*(tY«Sl”UY  )/SY”RC»V) 

H  1=  1 

IFCSU.LT. C)  Ml=-T 
S  L=  A 1 S  (  SU  ) 

SL=Sa*((Y-Sl-UY)/SY-RC«V) 

M2=l 

IFCSL.LT.C)  M 2  =  ” t 
S  L  =  A  °  S  C  SL  J 
TL=l./Cl. *P*SL) 

TL=  1-/C  1. *P*  SL  ) 

F  L= l”CA 1*TU  ♦  A  2*  HI**  2 ♦ A  3  *  T  U*  *3  ♦  A A*TU**4*  A5*TU*»5 )*E  XPC-SU**?) 
FL=l-(Ai*TL»A2*TL**?+A3»TL**3*A4*rL**4*A5*TL**5)*EXPC-SL**2) 
F=FXP(”  . 5*V  ** 2  )  •CMl*FU-M2*FL  ) 

FF=  l.C 

IFCCJ.EQ.  1T.0R.CJ.E0.2C1))  FE  =  0.E 

PF=PR*FF*F*CL 

1C  CONTINUE 

P»=SORT(?.)*PR/CA*SaRI(3.14I5S)) 

RETURN 

ENO 


A- 1 
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